Michael Ruggiero
# Function imports
import numpy as np
import pandas as pd
import networkx as nx
import osmnx as ox
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
ox.config(use_cache=True, log_console=True)
ox.__version__
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
The recent Merrimack Valley gas explosions reveled the need for cities to quickly deploy first responders while accounting for damaged roads. However, information during emergencies can be incomplete and properly dispatching police requires multiple information streams.
This project attempts to address this problem with the following workflow:
(1) Incorporate DOT api information about construction and accidents along flagged roads.
(2) Incorporate Twitter and Google information to estimate traffic patterns.
(3) Simulate a disaster zone in Medford that destroys a percentage of roads and creates emergency locations.
(4) Optimize dispatch routes for police to reach each disaster location.
# get a graph for some city
place = {'city' : 'Medford',
'state' : 'MA',
'country' : 'USA'}
medford = ox.graph_from_place(place, network_type='bike')
fig, ax = ox.plot_graph(medford)
# what sized area does our network cover in square meters?
medford_proj = ox.project_graph(medford)
nodes_proj = ox.graph_to_gdfs(medford_proj, edges=False)
graph_area_m = nodes_proj.unary_union.convex_hull.area
print(graph_area_m)
# show basic stats about the network
medford_stats = ox.basic_stats(medford_proj,
area=graph_area_m,
clean_intersects=True,
circuity_dist='euclidean')
# medford_stats = ox.extended_stats(medford_proj, ecc=True, bc=True, cc=True)
pd.DataFrame(medford_stats).set_index("streets_per_node_counts")
# Edge and Node projection
nodes_med, edges_med = ox.graph_to_gdfs(medford, nodes=True, edges=True)
nodes_med.head(2)
nodes_med.shape
#Functions built for this project
# %load graph_functions.py
import graph_functions as gf
medford_data = gf.node_roader(medford)
medford_data.keys()
edges_med, nodes_med, medford = medford_data["edges"], medford_data["nodes"], medford_data["graph"]
nodes_med.major_inter.value_counts()
#Long/Lat information
nodes_med[["x","y"]].describe()
nodes_med.head(2)
#A function to generate a random lat/long point
gf.random_point()
random_point_A = gf.random_point()
random_point_B = gf.random_point()
# get the nearest network node to each random point
orig_node = ox.get_nearest_node(medford, random_point_A)
dest_node = ox.get_nearest_node(medford, random_point_B)
# find the route between these nodes then plot it
route = nx.shortest_path(medford, orig_node, dest_node, weight='length')
fig, ax = ox.plot_graph_route(medford, route, node_size=0)
# Length of route in meters
nx.shortest_path_length(medford, orig_node, dest_node, weight='length')
# Absolute Distance between two nodes (as the crow flies)
ox.great_circle_vec(medford.node[orig_node]['y'], medford.node[orig_node]['x'],
medford.node[dest_node]['y'], medford.node[dest_node]['x'])
d_location = nodes_med[["y","x"]].sample(1).values[0]
print(d_location)
#returns a graphs object of a disater
disaster = gf.disaster_generator(nodes_med, location_point = d_location, radius = 2500)
nodes_dis, edges_dis = ox.graph_to_gdfs(disaster)
#Set Colors
ec = ['lightcoral' if i in disaster.edges() else 'green' for i in medford.edges()]
nc = ['red' if i in disaster.nodes() else 'blue' for i in medford.nodes()]
#Plot energency grid
ox.plot_graph(medford, node_size=15, node_color = nc ,edge_color=ec)
nx.degree_histogram(disaster)
medford_disaster = gf.road_kill(disaster,
nodes_med,
edges_med,
kill_percentage = .8 )
# get the boundary polygons for neighboring cities, save as shapefile, project to UTM, and plot
place_names = ['Medford, MA, USA',
'Arlington, MA, USA',
'Somerville, MA, USA',
'Malden, MA, USA',
# 'Melrose, MA, USA',
'Winchester, MA, USA',
'Stoneham, MA, USA',
'Everett, MA, USA',
]
medford_area = ox.gdf_from_places(place_names)
fig, ax = ox.plot_shape(medford_area)
# highlight one-way roads
oneway = ox.graph_from_place(place, network_type='drive')
ec = ['r' if data['oneway'] else 'b' for u, v, key, data in oneway.edges(keys=True, data=True)]
fig, ax = ox.plot_graph(oneway, node_size=0, edge_color=ec, edge_linewidth=1.5, edge_alpha=0.5)
#Neighboor hoods buffered
neigboors_buffered = ox.gdf_from_places(place_names, gdf_name='neighboors', buffer_dist=250)
fig, ax = ox.plot_shape(neigboors_buffered, alpha=0.7)
# DF builder of start/end nodes with long/lat and times
# Already done, commenting out to avoid clutter in notebook
router = gf.random_zone_picker(nodes_med, 10)
#router.to_csv("more_zones.csv")
router.head(3)
# network from address, including only nodes within 5.5km along the network from city hall
neighbor = ox.graph_from_address(address='85 George P. Hassett Drive,\
Medford, MA 02155',
distance=5500,
distance_type='network',
network_type='walk')
# you can project the network to UTM (zone calculated automatically)
neighbor_projected = ox.project_graph(neighbor)
# Edge and Node projection
nodes_area, edges_area = ox.graph_to_gdfs(neighbor_projected,
nodes=True,
edges=True)
nc = ['blue' if i in medford.nodes() else 'green' for i in neighbor.nodes()]
ox.plot_graph(neighbor,
node_size=15,
node_color = nc);
neighbor_data = gf.node_roader(neighbor)
black = neighbor_data["nodes"]
black[black.color == "black"].head(2)
neighbor_data.keys()
nodes_area, edges_area, neighborhood = neighbor_data["nodes"], neighbor_data["edges"], neighbor_data["graph"]
for i in neighbor_data['major_intersections']:
print("Number of major nodes: ",
len(neighbor_data['major_intersections'][i]),
"\t with degree ",
i,)
neighbor_data['major_intersections'].keys()
# find the route between these nodes then plot it
route = nx.shortest_path(neighbor, 3809869827, 65339393, weight='length')
fig, ax = ox.plot_graph_route(neighbor, route, node_size=0)
_ = gf.city_node(nodes_med,"medford",nodes_area)
nodes_area.head()
gmap, complete_node, node_occurances = gf.directionsToDataframe("10.Google_directions_api_all6k.csv")
gmap.head()
google_top = gf.traffic_node(nodes_area, node_occurances)
g_list = list(google_top[google_top.traffic_importance > 887].osmid.values)
google_top[google_top.traffic_importance > 887]
# Edge and Node projection
nodes_area, edges_area = ox.graph_to_gdfs(neighbor_projected,
nodes=True,
edges=True)
g_size = [25 if i in g_list else 0 for i in neighbor.nodes()]
ox.plot_graph(neighbor,
node_size=g_size,
node_color = "red", edge_alpha= .5);
import time
import pickle
from scipy.optimize import linear_sum_assignment
#Populate police station
"""
Medford Police:
100 Main St, Medford, MA 02155
42.415811 -71.110152
"""
police_node = gf.nodefinder(nodes_med,
lat = 42.415811,
lng = -71.110152)
def monte_simulator(p_e_samples, police_location, simulations, name = ""):
main_start = time.time()
police_node = gf.nodefinder(nodes_med,
lat = police_location[0],
lng = police_location[1])
police_station = [police_node] * p_e_samples
p_e_samples *= 2
monte_simulation = {}
#Iteration though simulation
for sim in range(0,simulations):
#Initialization
start = time.time()
print("Starting simulation: {} at time: {}".format(sim,round(start - main_start,2)))
monte_dis = {}
#Generates disaster and radius
disaster_location = nodes_med.sample(1)
monte_dis["dis_location"] = disaster_location
monte_dis["dis_rad"] = np.random.uniform(low = .4, high = .9) * 4000
monte_dis["disaster"] = gf.disaster_generator(nodes_med,
location_point = disaster_location[["y","x"]].values[0],
radius = monte_dis["dis_rad"],
plotter = 0)
monte_dis["nodes_effected"] = len(monte_dis["disaster"].nodes())
monte_dis["roads_effected"] = len(monte_dis["disaster"].edges())
#Randomly destroys roads in disaster radius
monte_dis["remain_per"] = np.random.uniform(low = .35, high = 1)
monte_dis["remaining"] = gf.road_kill(monte_dis["disaster"],
nodes_area, #These come from the neighborhood
edges_area,
kill_percentage = monte_dis["remain_per"],
plotter = 0)
#Make dictionaries for emergencies and police
e_list = np.random.choice(list(monte_dis["disaster"].nodes()),
replace=False,
size = p_e_samples)
p_list = list(nodes_med.sample(int(p_e_samples / 2)).osmid.values)
p_list.extend(police_station)
monte_dis["emergency"] = {}
monte_dis["patrol"] = {}
for i in range(0,p_e_samples):
monte_dis["patrol"]["officer_{}".format(i)] = p_list[i]
monte_dis["emergency"]["emergency_{}".format(i)] = e_list[i]
rows = list(monte_dis["patrol"].keys())
columns = list(monte_dis["emergency"].keys())
#Distance Matrix
d_matrix = pd.DataFrame(columns = columns,
index = rows)
for i in range(len(rows)):
for j in range(len(columns)):
try:
path = nx.shortest_path_length(monte_dis["remaining"],
monte_dis["patrol"][rows[i]],
monte_dis["emergency"][columns[j]],
weight='length')
except:
path = 10_000
d_matrix.loc[rows[i], columns[j]] = path
monte_dis["distance_matrix"] = d_matrix.astype(float)
#Use Hungarian algorithm for linear sum assignment
lsa = linear_sum_assignment(monte_dis["distance_matrix"].to_numpy())
adj_matrix = np.matrix(np.zeros((10, 10)))
#optimized distances
od_list = []
for i in range(len(lsa[0])):
od = monte_dis["distance_matrix"].iloc[lsa[0][i],lsa[1][i]]
adj_matrix[lsa[0][i],lsa[1][i]] = od
od_list.append(od)
dispatch = {"officer_{}".format(lsa[0][i]):["emergency_{}".format(lsa[1][i]),
od_list[i]] for i in range(len(lsa[0]))}
monte_dis["ad_matrix"] = adj_matrix
monte_dis["dispatch"] = dispatch
monte_dis["dispatch_distance"] = sum(od_list)
monte_simulation[sim] = monte_dis
if sim % 200 == 0:
print("pickle dump")
pickle.dump(monte_simulation, open( "monte_simulation_{}_{}.p".format(name,sim), "wb" ))
end = time.time()
print("Ending simulation: {}. Total {} seconds for run\n".format(sim,round(end - start,2)))
return monte_simulation
monte_simulation = monte_simulator(5,[42.415811,-71.110152], 5, "demo")
monte_simulation.keys()
monte_simulation[0].keys()
monte_simulation[0]['dispatch']
monte_simulation[2]['dispatch']
monte_simulation[2]["ad_matrix"]
size_list = range(len(monte_simulation))
df = pd.concat([monte_simulation[i]["dis_location"].copy() for i in size_list]).reset_index(drop=True)
extra_data = ['dis_rad', 'nodes_effected', 'roads_effected', 'remain_per', 'dispatch_distance']
for i in range(len(monte_simulation)):
for j in extra_data:
df.loc[i, j] = monte_simulation[i][j]
df.columns
X = df[['osmid','x', 'y','dis_rad', 'nodes_effected', 'roads_effected',
'remain_per', 'dispatch_distance']]
X.head()